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Abstract 

Many viruses rely on the self-assembly of their capsids to protect and transport their genomic 
material. For many viral systems, in particular for human viruses like hepatitis B, adeno or human 
immunodeficiency virus, that lead to persistent infections, capsomeres are continuously produced 
in the cytoplasm of the host cell while completed capsids exit the cell for a new round of infection. 
Here we use coarse-grained Brownian dynamics simulations of a generic patchy particle model to 
elucidate the role of the dynamic supply of capsomeres for the reversible self-assembly of empty 
T1 icosahedral virus capsids. We find that for high rates of capsomere influx only a narrow range 
of bond strengths exists for which a steady state of continuous capsid production is possible. For 
bond strengths smaller and larger than this optimal value, the reaction volume becomes crowded 
by small and large intermediates, respectively. For lower rates of capsomere influx a broader range 
of bond strengths exists for which a steady state of continuous capsid production is established, 
although now the production rate of capsids is smaller. Thus our simulations suggest that the 
importance of an optimal bond strength for viral capsid assembly typical for in vitro conditions 
can be reduced by the dynamic influx of capsomeres in a cellular environment. 


* These authors contributed equally. 
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I. INTRODUCTION 


Viruses are experts on the inner working of cells and their investigation has contributed 
strongly to our understanding of the physical principles at work in biological systems. In par¬ 
ticular, the study of viruses has demonstrated the amazing power of biological self-assembly 
in an experimentally and theoretically accessible system. In their physiological context, 
viruses rely on the molecular machinery of their host to reproduce both their genomic mate¬ 
rial and the protein capsid usually encapsulating it jl, 2], Typically the capsid is assembled 
from many copies of only a few different capsid proteins and shows icosahedral or helical 


symmetry 


-7]. The elementary assembly blocks for the assembly of a capsid are termed 


capsomeres. They can either consist of single capsid proteins or of preassembled sets of 
capsid proteins. For many viruses the formation of the capsid can be reproduced in vitro 
[<3, 9]. The dynamics of in vitro assembly has been analyzed using light and small-angle 
X-ray scattering techniques 10H13I]. While capsid assembly of viruses with single-stranded 
genomic material often requires the presence of the genomic material as an ’electrostatic 
glue’ [l, 14], this is typically not possible for double-stranded genomic material due to its 
larger bending stiffness. These viruses typically assemble their capsid without the genomic 


Im¬ 


material, which is then inserted into the capsid by a motor {l 

Despite the plethora of known capsid structures [l(3], the dynamic assembly process of 
the virus shell is still far from being fully understood. As experimental possibilities for de¬ 
tailed monitoring of the assembly dynamics are limited, modeling can significantly help to 
increase our understanding of the mechanisms that govern the assembly process. In the past 


various techniques ranging from coarse-grained molecular c 


simulations 
proaches 


17 


29 


251 ] through Monte Carlo simulations 


331] to thermodynamic descriptions 


10 


34 


26 


ynamics or Brownian dynamics 
and discrete stochastic ap- 


361 ] have been used to elucidate dif¬ 
ferent aspects of the assembly process. For example, the influence of capsomere shape [ 19 I. 24] 
or the emergence of polymorphic structures [37H39] have been investigated with such theo¬ 
retical approaches. The different techniques used to study virus capsid assembly have been 
recently reviewed by Hagan 40]. 

Several studies have shown that the successful assembly of complete capsids starting 
from a fixed number of assembly subunits requires intermediate (or even optimal) bond 


strengths 
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42], At low bond strengths (or high temperature) 
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the formation of large clusters is suppressed. At high bond strengths (or low temperature), 
by contrast, the simultaneous formation of stable assembly intermediates, which cannot 
recombine into complete capsids, or the formation of large, misassembled structures with 
non-native binding interactions can prevent the formation of complete capsids. In this case 
the system becomes kinetically trapped. 

While in vitro studies and computer simulations usually work with a finite initial number 
of capsomeres that are increasingly used up during the capsid assembly process, in the 
physiological context of the cell the capsomeres are produced in a continuous fashion. In 
a recent stochastic simulation study, it has been shown for the case of genome-stabilized 
capsids (as for example for the bacteriophage MS2) that this ’protein ramp’ can make virus 
assembly very robust against kinetic trapping 43]. While bacteriophages typically kill the 


host cell to start a new round of infection, many animal and plant viruses tend to follow 
strategies that keep the host cell alive at least for a certain period of time. This is especially 
true for human viruses such as hepatitis B virus (HBV), human adeno virus (HAdV) or 
human immunodeficiency virus (HIV), which lead to persistent infections. Recently the 
role of dynamic protein supply for viral capsid assembly has been studied from a systems 
level perspective using a kinetic gene expression model with exponential protein production 


and a master equation for capsid assembly 


44], However, the effect of continuous protein 


production has not been studied yet in a spatial model. For empty capsid assembly, it has 
been shown earlier with Brownian dynamics that steady states exist in which the removal of 


large clusters is balanced by the reinsertion of the corresponding monomers 


42], However, 


the effect of the capsomere supply rate on these steady states has not been studied yet. 

Inspired by the notion of continuous virus production, we use coarse-grained Brownian 
dynamics simulations to investigate the assembly of empty T1 capsids in the presence of 
a dynamic capsomere supply. In order to avoid crowding of the reaction volume by large 
clusters and motivated by the exit of completed virions from the cell by budding or exocyto- 


sis 


45 


48], capsids are removed from the simulation in the moment that they are completed. 


Our simulations suggest that there exists a certain range of bond strengths in which the in¬ 
flux of new capsomeres is balanced by capsid removal. This steady-state region is surrounded 
by parameter regions in which successful capsid assembly is prohibited by crowding of the 
reaction volume; depending on influx rate and bond strength, crowding is observed either 
with small or large intermediates. Our main result is that the favorable region for continuous 
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FIG. 1. Schematic illustration of our simulations. Capsomeres are randomly inserted into the 
simulation volume with a rate ki while completed capsids are removed immediately (including 
the intermediates they might contain inside). This simple model mimics the situation in human 
cells with persistent infections in which capsomeres are continuously produced by translation and 
completed capsids leave by budding or exocytosis. 

virus production becomes larger for lower influx rates. Our work identifies essential limits 
of viral self-assembly in a dynamic context and suggests that bond strength has to be less 
fine-tuned in a cellular context than in in vitro experiments. 


II. METHODS 


To investigate the effect of a continuous influx of capsomeres on virus assembly, we use 
an efficient coarse-grained Brownian dynamics approach, which has previously been used to 


study the effect of reactivity switching during the assembly process [25]. Here we consider a 


T1 capsid which is composed of 60 identical capsomeres ( 7 ]. Each capsomere is described by 
a hard sphere which is equipped with spherical patches (see Fig.[U). The spatial arrangement 


of the patches re 


Berger et al. 


17 


iects the capsid geometry according to the local rules scheme developed by 
49]. All assembly intermediates formed during the assembly process are 


treated as rigid objects and are propagated according to their translational and rotational 


diffusive properties 


50 


5l|, which are evaluated on-the-fly upon their formation 
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If an overlap between two complementary patches is realized by diffusion, a bond is formed 
with the probability P reac t = k a At <C 1. This bond can be established either between two 
unconnected clusters (inter-bond) or in an already connected cluster (intra-bond). Otherwise 
we only consider one type of bond to keep the number of parameters small. We assume that 
bond formation is achieved by strong local forces such that this process is very fast on 
the time scale of our simulations. Therefore upon formation of an inter-bond, the clusters 
instantaneously assume the correct relative position and orientation for the assembly of the 
capsid unless the necessary reorientation results in a steric overlap either between the two 
merging clusters or with other clusters. We note that our approach based on patchy particles 
and local rules does not allow us to study the formation of aberrant cluster with non-native 


interactions 


20 


22 


42], Because we do not consider any forces, our approach also does 


not allow us to investigate strained capsids, as possible in molecular dynamics or Brownian 
dynamics simulations with potentials. 

In order to study reversible dynamics, every existing bond can also dissociate with the 
probability P dissoc = k d At < 1. If bond dissociation results in two unconnected clusters, 
they are positioned relative to each other according to a computational scheme which ensures 
detailed balance in order to prevent additional, non-physical driving forces for the self- 
assembly [53]. 

Intra-bond formation in an already connected cluster leads to an additional stabilization 
of the cluster as closed loops are formed, in which every capsomere is connected to at least 
two neighboring capsomeres. For such a loop structure to break apart it is necessary that 
all bonds have to be in the open state simultaneously. The formation/dissociation of an 
intra-bond does not affect the structure of the complex. The energy gain by bond formation 
is related to the microscopy rates by E = —ln(/c a //c d ), where k b is the Boltzmann 


constant 


53]. 


To study the effect of a continuous supply of capsomeres on the assembly dynamics, we 
introduce the influx rate k n . We place a new capsomere in the simulation volume in each time 
step with the probability p in = kiAt -C 1. Position and orientation of the new capsomere 
are randomly chosen with the constraint that no steric overlap with existing clusters is 
created. Inspired by virus exit strategies like budding or exocytosis, which in principle 
can lead to a steady state of virus production, complete capsids are removed immediately 
upon formation together with all intermediates inside the capsid (compare Fig. [[]). This 
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removal rule is based on the assumption that complete capsids are much more stable than 
partially assembled capsids and that cellular mechanisms exist that are exploited by the 
virus to leave the cell. We note that our model does neither incorporate any details of 
the production mechanism nor any details of the exit mechanism, but is kept as simple as 
possible in order to investigate the underlying physical principles of the assembly process in 
such a dynamic setup. For the same reason we focus on the assembly of T1 virus capsids 
and assume all bonds to be identical. 

For the following it is helpful to introduce the concept of dimensionless bond strength 
k s , which is defined by the ratio of the microscopic association rate k a and the microscopic 
dissociation rate k d : 


k s = K/k d . 


( 1 ) 


Note that bond strength is similar to, but different from the equilibrium association constant 
K eq for a bimolecular reaction, because it is defined by the ratio of microscopic rates (with 
the physical dimension 1/s) rather than by the ratio of a macroscopic association rate 
constant k on (with physical dimension l/(sM)) and a macroscopic dissociation rate k a s 
(with physical dimension 1/s). For the reaction between two clusters (without any closed 
loops) K eq = V*k s is related to k s by the encounter volume V* (with physical dimension 
m 3 ) which is defined by all two-particle configurations of the two clusters with an overlap of 


complementary patches 


53]. For the case of virus assembly investigated here, k s is not only 


a measure for the strength of inter cluster bonds, but also for the stability of closed loops. 

In order to avoid the need to exhaustively scan parameter space and motivated by previous 
results quantifying the success of assembly as a function of k a and k d |25j, we use a linear 
relation between dissociation and association rate k d = mk a + c with c = 0.0111ns -1 and 
m = —0.0011 to explore the parameter space ranging from very strong to very weak bond 
strength. When varying k s from 10 2 to 10 5 we explore association and dissociation rates in 
the range of k a £ [Ins -1 ,10ns -1 ] and k d £ [1 x 10 -4 ns -1 , 0.01 ns -1 ], respectively. 

All simulations have been performed at a time resolution of At = 0.01ns using periodic 
boundary conditions. Capsomeres are modeled as hard spheres of radius R s ter i C = 1 nm. 
Each capsomere is equipped with three distinct spherical patches reflecting the geometry 
of the T1 capsid. Each patch has a radius of r pat ch = 0.3 nm with the center of the patch 
being located on the surface of the hard sphere. A bond can only be established between 
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complementary patches according to the local rules. The diffusive properties of all interme¬ 
diates are represented by their mobility matrices which are evaluated at room temperature 
T = 293 K using the viscosity of aqueous medium rj — 1 mPa s. 

Given the typical time and length scales Iq — 1 nm and to = 1 ns of our simulations, we 
define the dimensionless length parameter A = l/l 0 and the dimensionless time parameter 
r = t/to in order to simplify the notation. Furthermore, we introduce the dimensionless 
box volume A = VboxAo, the dimensionless particle concentration p = IV/A, where N is 
the number of capsomeres in the simulation volume, and the normalized influx rate Hi = 
fcjtol0 6 /A. Hi can be understood as the rate of concentration increase due to the influx 
of capsomeres. Using the normalized influx rate Hi instead of ki allows us to compare the 
assembly process for different sizes of the simulation volume. Our simulation volume has a 
typical linear extension of 30 to 50 nanometers, leading to reasonable computing times for 
complete capsid assembly. 

Considering particles which are equipped with only one spherical patch we can relate the 
values of k s used here to K cq and estimate the speed of the reaction dynamics for this case. 
In contrast to the capsomeres these particles can only form dimers. The encounter volume 
for this reaction is V * ~ O.llnm 3 (compare reference 53] for details on the calculation) and 
the equilibrium association constants K eq range from 6.6M -1 to 6.6 x 10 3 M -1 for the values 
of bond strength used here (10 2 < k s < 10 5 ). Thus our simulations proceed at a relatively 
high concentration in the mM-range. 


In order to estimate the influence of the pate 


algorithm developed by Zhou and coworkers 


i size for the speed of the reaction, we use an 


54 


55[ with which the diffusive association rate 


constant kp can be estimated based on the survival probability of two clusters starting in 
an encounter 53]. For the dimerization we estimate kp ~ 6.25 x 10 8 M _1 s _1 . The diffusive 
dissociation rate constant then follows as kn,b — kjj/V *. Depending on the values of the 
microscopic reaction rate k a the macroscopic association rate constant k + = k D k a /(k a + 
kn,b) is in the range of 6.0 x 10 7 M -1 s -1 to 3.21 x 10 8 M _1 s _1 . For lower values of k + 
the reaction can be considered as reaction-limited {k a < ko,b) while for higher values of 


k + the assembly is equally influenced by reaction and diffusion {k a ~ ko,b) [56]. Although 


bimolecular reactions in this range of macroscopic association rate constants k + are typically 
considered diffusion-limited 57], the relatively large patch sizes used here allow for rapid 


formation of the diffusive encounter so that the reactions in our simulations are at least 
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partly reaction-limited according to the classification scheme by Eigen 56]. The macroscopic 
dissociation rate k- = kD,bkd/{k a + ko,b) ranges from 4.8 x 10 4 s _1 to 9.0 x 10 6 s _1 , which is 
approximately the same range as for the microscopic dissociation rate. 

It is instructive to compare the reaction rate constants and equilibrium constants for 
dimerization used in our simulations to the rates inferred by Xie at al. 58| with a non-spatial 
model from published light-scattering data of the in vitro assembly of hepatitis B virus 
(ffBV) virus 10], cowpea chlorotic mottle virus (CCMV) [ll] and human papillomavirus 
(ffPV) 12], We first note that Xie et al. consider closed loops as being infinitely stable 
which allows for capsid assembly at lower concentrations than in our case. Nevertheless 
we find that for ffBV and CCMV, our equilibrium association constants are well within 
the range reported by Xie at ah, while our association rate constants are larger by two 
orders of magnitude. In this context, we note that the use of high concentrations and 


enhanced assernb 
assembly 


20 ], 


y dynamics is a common limitation of particle-based simulations of capsid 


23|, 


25 


59 


60|. 


In our setup, concentrations are not fixed, but arise from the dynamic influx. In practice, 
capsomere and capsid production rates can vary widely for different viruses and different 
host systems. For influenza virus, for example, a production rate of 10 4 virions over 10 
hours has been reported after start of viral protein translation 61]. For HIV, in contrast, 


only 800 virions are produced over 8 h after start of viral protein translation 


62], Similar 


variations also exist for the transcription and translation rates. As we will see below, our 
capsid production rate is of the order of 10 8 virions per hour in a small reaction volume 
because we consider a small generic virus with only 60 identical components, relatively high 
influx rates and accelerated dynamics. 


III. RESULTS 

First we qualitatively characterize typical responses of our simulation setup at different 
bond strengths. For this purpose the time course of the total number of capsomeres placed 
inside and removed from the simulation volume is shown in Fig. [2] for two different values 
of k s and a normalized influx rate of k* = 2.593 (a movie for the inital stages of capsid 
assembly is provided as supplementary movie SI). The number of capsomeres placed inside 
the simulation volume (black, solid line) shows the expected linear increase defined by the 





















FIG. 2. Evolution of the total number of capsomeres placed inside and removed from the simulation 
volume. The black, solid line depicts the number of capsomeres placed inside the simulation volume 
with Ki = 2.593. The red, dashed lines show individual trajectories of the number of capsomeres 
which are removed from the simulation volume due to capsid completion for a bond strength of 
k s = 1179, resulting in a steady state. The blue, dotted lines show individual trajectories of 
removed capsomeres for a bond strength of k s = 8483 at which crowding occurs. 

influx rate. For the low bond strength of k s = 1179 (red, dashed lines) two individual 
trajectories are shown. For these trajectories we see that after an initial lag phase without 
capsid completion (the length of the lag phase depends strongly on the rate of capsomere 
influx and the established concentration in the simulation volume (0.5 —2.0 x 10 4 ns)) a steady 
state with constant capsomere concentration is established in which capsomere influx and 
capsid removal balance each other (a movie for the steady state is provided as supplementary 
movie S2). 

For the higher bond strength of k s = 8483 (blue, dotted lines) the trajectories behave very 
differently. After the initial lag phase the rate of capsomere removal due to capsid completion 
is almost compensating the capsomere influx for a certain period of time until the rate of 
capsid completion drastically slows down. The time point at which this slow-down in capsid 
production is observed strongly varies between different trajectories. Once capsid completion 
has started to slow down, the concentration in the simulation volume quickly increases due 
to the influx of further capsomeres and the simulation volume eventually becomes crowded. 
If the volume fraction of the simulation which is occupied by capsomeres is too large, further 
assembly is prohibited. For the following we therefore introduce a crowding-threshold. This 
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FIG. 3. Phase diagrams with the different assembly regimes for different simulation volumes, (a) 
and (b) show the averaged capsomere concentration p as a function of influx rate Ki and bond 
strength k s for a simulation volume of Ai = 27000 and A 2 = 42875, respectively, (c) and (d) show 
the fraction of trajectories which became crowded during the simulation time as a function of Ki 
and k s for a simulation volume of Ai and A 2 , respectively. Each data point is obtained from 16 
independent trajectories using a simulation time of r s i m = 10 6 . 

threshold is reached when 35% of the simulation volume is occupied by capsomeres and 
trajectories reaching this threshold are aborted. The chosen threshold is much larger than the 
highest steady-state concentration established in our simulations and a trajectory reaching 
this high concentration will certainly lead to a full stop of assembly. 

We now systematically investigate the assembly process as a function of bond strength 
k s and normalized influx rate Ki. In our dynamic setup the rate of capsid production is pre¬ 
scribed by the influx of capsomeres and hence cannot be used as a measure for the quality of 
the assembly process. Instead we measure the ability of the assembly process to compensate 
the influx of new capsomeres. As can be seen for the case of a steady state in Fig. [5] the con¬ 
centration in the simulation volume linearly increases until a steady concentration is reached 
at which the assembly of full capsids compensates the influx of new monomers. Thus, for a 
driven system with a continuous influx of capsomeres the concentration established in the 
simulation volume for a given Ki is a measure for the efficiency of the assembly process: the 
lower the established concentration, the higher the efficiency (ability of the assembly process 
to compensate the influx). 

In Fig. [3] (a) and (b) the average capsomere concentration p during our simulations 
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is shown for two different simulation volumes, respectively. Here only those parameter 
combinations are shown for which a steady state is established in all trajectories during the 
simulation time. If no steady state is established (X), the average concentration does not 
provide a valid measure for the efficiency of the assembly process as the concentration in 
trajectories showing crowding rapidly increases until they are finally aborted. Hence for 
parameter combinations where at least one trajectory shows crowding we instead quantify 
the ability of the assembly process to compensate the influx of new capsomeres by the 
fraction of crowded trajectories <~p\ the higher the crowding fraction is the smaller the ability 
of the assembly process to compensate the influx of new capsomeres. This is shown in 
Fig. [3(c) and (d) for the two different simulation volumes, respectively. 

In general, we see from Fig. [3] (a) and (b) that for every value of k* used here the 
lowest concentration is found at intermediate bond strengths of k s ~ 10 3 . This is similar to 
assembly under static conditions where intermediate bond strengths show the largest yield 
of full capsids. Moreover, a minimum bond strength of around 2 x 10 2 is necessary for any 
assembly to occur. Below this threshold in bond strength the simulation volume becomes 
crowded almost independent of the normalized influx rate 

For high normalized influx rates we see that only a narrow range of intermediate bond 
strengths exists for which a steady state is established in our simulations in all trajectories, 
and the steady-state concentration even at optimal k s is very high. Here the simulation 
volume does not only become crowded for very small k s , but a distinct crowding regime 
exists at high k s . Thus, in the case of a high forcing of the system due to a rapid influx 
of capsomeres (corresponding to a fast production of capsids) establishing a steady state 
requires an optimal bond strength. 

By decreasing the normalized influx rate of capsomeres (lower forcing of the system) 
the average concentration established in the simulation volume also decreases (Fig. [3] (a) 
and (b)) as there is more time for the assembly to proceed in between the addition of 
new capsomeres. Interestingly, the crowding regime at high bond strengths shrinks with 
decreasing Ki (Fig. [31(c) and (d)), and for /q : = 1.111 a broad range of k s (also extending to 
very high k s ) exists for which a steady state of capsid production is observed. This is a clear 
difference between the dynamically driven system analyzed here and static systems with a 
fixed capsomere concentration. While in the static case kinetic trapping would prevent any 
capsid assembly at high values of k s , the continuous influx of new capsomeres reduces the 
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requirement for an optimal bond strength and successful assembly is possible even at very 
high k s . 

As Ki has been normalized by the simulation volume and can be understood as the rate 
of increase in concentration due to the influx of capsomeres, one expects our findings to 
be independent of the system size. Comparing the average concentration p in the smaller 
simulation volume (Fig. [31(a)) and in the larger volume (Fig. [31(b)) we indeed see that the 
average concentrations are almost identical in both cases if the rate of capsomere influx 
is appropriately scaled. This allows for a volume independent comparison of the effect of 
a dynamic influx of capsomeres on the assembly process. On the contrary, the crowding 
tendency especially at high bond strengths and high values of K t is reduced in the larger 
simulation volume (Fig. Eli) when compared to the smaller simulation volume (Fig.[3b). Thus 
the crowding tendency depends on the system size in a non-trivial manner and reflects the 
stochastic nature of the crowding process with fluctuations being reduced in larger systems 
as will be discussed in detail below. 

In order to further characterize the different regimes, individual trajectories of the time 
evolution of the capsomere concentration p for three different bond strengths with = 2.592 
and A = 27000 are shown in Fig. [4] (a)-(c) (the lower row shows corresponding simulation 
snapshots). Fig. Oja) shows the evolution of p of independent trajectories for intermediate 
bond strength k s = 1179. Here no crowding of the simulation volume occurs and a steady 
state is established in all simulations. After an initial phase with an increase in p a constant 
concentration is established in all simulations and individual trajectories stochastically fluc¬ 
tuate around this constant concentration. Thus, the average concentration in this case can 
indeed be used to characterize the system. 

In Fig. SKb) the time evolution of p for low bond strength k s = 164 is shown. In this case 
the bond strength is below the minimum k s needed for assembly of larger clusters, and the 
concentration in all trajectories linearly increases with simulation time leading to an almost 
deterministic abortion of the simulations due to crowding. 

In Fig. SKc) the evolution of concentration is shown for high bond strength k s = 8483. 
This case corresponds to crowding of the simulation volume at high bond strengths k s and 
high normalized influx rate Kj. Compared to the almost deterministic crowding process at 
very low bond strengths (Fig. ED(b)), we see that the crowding characteristic for this case is 
fundamentally different as now crowding of the simulation volume occurs stochastically. In 
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FIG. 4. Time evolution of the capsomere concentration for different assembly regimes («j = 2.592, 
A = 27000). For each of the different regimes one representative snapshot is shown in the lower 
row, with size-dependent color coding. Simulation data for (a) assembly without crowding at 
intermediate bond strength k s = 1179, (b) assembly with crowding at low bond strength k s = 164, 
and (c) assembly with crowding at high bond strength k s = 8483. The dashed black line represents 
the crowding limit at which trajectories are aborted. 

all trajectories shown in Fig. a very high, quasi-const ant concentration is established for 
a certain period of time until the system stochastically reaches an unfavorable configuration. 
Once such a configuration is reached, the formation of further capsids is hindered. This leads 
to an increase in concentration, which further slows down capsid production presumably due 
to steric collisions during the reaction process, and the system then quickly becomes crowded 
with large capsid intermediates due to the addition of new capsomeres. 

A similar crowding characteristic is observed at high K t for bond strengths (k s m 268.0) 
which are slightly above the threshold in k s necessary for any capsid assembly. In this case 
again capsid completion can balance capsomere influx for a certain period of time until the 
simulation volume quickly becomes crowded albeit in this case with a high fraction of small 
clusters. 
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FIG. 5. Relative population of different cluster sizes pk for k ? ; = 2.593 and different values of k s 
as a function of the cluster size. The red, solid lines correspond to bond strengths without any 
crowding of the simulation volume (439 < k s < 1931). The blue, dashed and dotted lines and the 
green, dashed lines correspond to bond strengths where all trajectories showed crowding at high 
(k s > 8483) or at low bond strengths (k s = 100 and k s = 164), respectively. 

In general, stochastic crowding of the simulation volume requires a very high concen¬ 
tration of capsomeres. In this case the realization of an unfavorable configuration hinders 
capsid assembly for a certain time and the influx of new capsomeres leads to a further in¬ 
crease in concentration which in turn hinders further assembly due to steric collisions. Thus, 
once a certain concentration is surpassed in our simulations the simulation volume inevitably 
becomes crowded as the assembly process cannot compensate the influx of new capsomeres. 
The observation that the crowding process is triggered stochastically by the realization of 
an unfavorable configuration of the system agrees well with the previous observation that 
the stochastic crowding tendency decreases with system size (compare Fig. [3] (c) and (d)). 
For a larger system the relative fluctuations in concentration decrease and the probability 
that an unfavorable configuration of the whole system is realized is reduced. 

In order to further analyze the mechanisms leading to crowding we compare the relative 
population of cluster sizes pk = fk x k/N for different cases. Here fk is the number of 
clusters of size k and N the total number of capsomeres in our simulation volume. Thus pk 
is the probability that an arbitrarily chosen capsomere is part of a k-sized intermediate. In 
Fig. [5] pk is shown for different bond strengths and a normalized influx rate of Kj = 2.593 
on a logarithmic scale. Here pk has been first averaged over the whole simulation time of a 
trajectory and subsequently over 16 independent trajectories. 
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For the intermediate bond strengths for which no crowding is observed (red, solid lines) 
small and large cluster sizes dominate. This is in agreement with previous observations char¬ 
acterizing successful assembly from a fixed concentration of capsomeres 19, [hi] and indicates 
that in this regime successful assembly proceeds by the addition of small clusters to only 
a few larger, stable assembly intermediates. For crowded runs at very low bond strengths 
(green, dashed lines) only small cluster sizes are populated and no capsid completion is ob¬ 
served. This explains the quasi-deterministic increase in concentration shown in Fig. (4](b). 
For crowded runs at very high bond strengths (blue, dashed and dotted lines) the relative 
population at large cluster sizes strongly resembles the relative population without crowd¬ 
ing. In contrast the population of small cluster sizes is strongly reduced for the crowded 
runs while the population of intermediate cluster sizes (k m 10) is increased. This shows 
that in this case new capsomeres are quickly absorbed by existing clusters and the system 
becomes eventually crowded with intermediate and large capsid fragments. Interestingly 
the quick depletion of single capsomeres has previously been identified as characteristic for 


kinetic trapping during capsid assembly 25 


35 


36 


41| . This shows that kinetic trapping 


and box crowding at high bond strengths are strongly connected phenomena. 


As we have seen, crowding effects are essential to understand capsid assembly, but until 
now we only have considered self-crowding by viral components. In the cell, crowding will be 
established also by other crowders and therefore lower (more physiological) concentrations 
of viral components are expected to be sufficient to result in similar effects as described 
here. In order to test the validity of our findings in the presence of additional, non-specific 
macromolecular crowders, we have used our dynamic simulation setup with non-reactive 
crowders. To this end we place a total of 600 spherical crowders inside a simulations volume 
of size Ai = 27000. These crowders have the same radius as the capsomeres, however, they 
do not participate in any reactions (no patches). In the simulation volume we now have a 
total concentration of ptotai = p + Powder with the concentration of crowders p cr owder being 
kept constant throughout the simulations. This implies that any crowder found inside a 
complete capsid is placed back into the simulation volume. The upper (red) histogram in 
Fig. O shows the dependence of the average concentration p observed during our simulations 
on the bond strength k s for a normalized influx rate of K; = 0.185. Again p is only shown for 
those values of k s for which a steady-state was established in all trajectories during the sim¬ 
ulation time. The lower (blue) histogram depicts the fraction of crowded runs <p for different 
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FIG. 6. Averaged capsomere concentration p and fraction of crowded trajectories i^asa function of 
bond strength k s for a normalized capsomere influx rate of m = 0.185 when placing 600 additional, 
non-specific crowders inside a simulation volume of size Ai = 27000. Each data point is obtained 
from 16 independent trajectories using a total simulation time of r s i m = 6 x 10 5 . On the right a 
representative simulation snapshot is shown for k s = 1179 and with size-dependent color coding of 
the capsid fragments. The additional crowders are colored in gray. 

bond strengths. In addition a representative snapshot of the assembly system including the 
additional crowding agents is shown in Fig. [6j Comparing the dependence of p on the bond 
strength k s , we see the same qualitative behavior with additional crowders (Fig. [6]) and 
without crowders (Fig. [3]). In both cases assembly is most efficient at intermediate bond 
strengths. Furthermore in Fig. Owe again observe two regions (at very low and very high 
bond strength) in which the simulation volume becomes too crowded for a steady state to 
be established in all trajectories. This is similar to the two self-crowding regions observed in 
Fig. [3] at high normalized rates of capsomere influx. However, the normalized influx rate of 
capsomeres used in Fig. [6] is smaller than those used in Fig. [3j This suggests that, although 
the qualitative dependence remains the same, it is more difficult to establish a steady-state 
virus production in already pre-crowded environments and lower rates of influx are required 
in this case. 

After having analyzed the different mechanisms which prevent a steady state from being 
established, we now focus on small normalized influx rates. As can be seen in Fig. [3](a) and 
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FIG. 7. Comparison of capsid yield for a fixed initial capsomere concentration and for a continuous 
influx of capsomeres. Here the yield of complete capsids in a simulation volume of A = 125000 
is shown as a function of the bond strength k s after a simulation time of r s i m = 10'. The blue 
histogram shows the yield when initially placing 1000 capsomeres in the simulation volume. The 
red histograms on the other hand show the yield of capsids when gradually increasing the number 
of capsomeres from zero with a rate of K{ = 0.24 (light red) and m = 0.08 (dark red), respectively, 
until a total of 1000 capsomeres is reached. The yield has been averaged over 16 different runs. 

(b) a steady state with continuous capsid production is established even for very high bond 
strengths if the normalized influx rate is small enough. This suggests that a gradual influx 
of capsomeres can prevent or at least reduce kinetic trapping. 

In order to verify that dynamic capsomere influx indeed reduces the requirement of an 
optimal bond strength we compare the yield of capsids for two different setups: a static 
setup initially starting with No = 1000 capsomeres and a dynamic setup in which a total 
of 1000 capsomeres is placed in the simulation volume with a certain rate. In both cases 
complete capsids are considered to be stable and are taken out of the simulation volume. In 
Fig. [7] the yield of full capsids after a simulation time of r sim = 10' is shown as a function 
of the bond strength k s . The blue histogram shows the yield in the static case while the 
red histograms show the yield when gradually placing new monomers into the simulation 
volume with a normalized influx rate of k* = 0.24 (light red) and k, = 0.08 (dark red). 

At low bond strengths ( k s < 720) the yield in the static case is slightly higher than 
in the dynamic case. Thus, for small bond strengths it is beneficial if a larger number of 
capsomeres is available throughout the whole simulation time. In the case of optimal bond 
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strengths (k s = 1179 — 1931) no difference between the two setups is observed as assembly 
at these bond strengths proceeds quickly, and almost all simulations show the maximum 
yield of 16 capsids within the simulation time. 

When further increasing the bond strength we see, however, that the yield in the static 
case quickly drops and no complete capsids are observed above k s = 8483. In this case 
kinetic trapping completely prevents the formation of capsids. For the dynamic case, in 
contrast, we still observe considerable yield of capsids above k s = 8483 with the yield being 
higher for the lower normalized influx rate used. This indeed shows that a dynamic setup 
with a gradual supply of capsomeres reduces the selectivity for an optimal bond strength 
and makes the assembly process more robust and less vulnerable against kinetic trapping. 


IV. CONCLUSION 


The assembly of the viral protein shell, the capsid, from elementary assembly units (cap¬ 
someres) is a key step during the replication of most viruses. The assembly process needs to 
be sufficiently robust to guarantee the successful formation of the capsid in the dynamic en¬ 
vironment of the host cell. While for test tube experiments on capsid assembly the material 
available for the assembly process remains constant, in the cellular environment the elemen¬ 
tary building blocks for the assembly process are continuously produced by the biomolecular 


machinery of the cell 


43] 


Here we have investigated the role of a dynamic supply of capsomeres and the removal 
of complete capsids for the assembly of empty T1 virus capsids by using a minimal spatial 
model based on coarse-grained Brownian dynamics simulations. It has been shown earlier 
that such a setup can result in a steady state with capsomere influx being balanced by 
capsid completion 42]. Our simulations reveal that for very high rates of capsomere influx 


the assembly process is only able to compensate the influx of new capsomeres in a narrow 
range of intermediate bond strengths while outside this range the simulation volume be¬ 
comes crowded. At lower bond strengths the formation of larger clusters is prohibited and 
the simulation volume becomes crowded almost deterministically with small clusters. At 
higher bond strengths, in contrast, the simulation volume becomes crowded with large, in¬ 
compatible assembly intermediates. This crowding process is triggered stochastically by the 
formation of an unfavorable configuration of the system. While crowding of the simulation 
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volume at very low bond strengths is nearly independent of the influx rate of capsomeres, the 
crowding regime at high bond strengths vanishes for a slower influx of capsomeres. Thus, 
for smaller rates of capsomere influx a steady state with continuous capsid production is 


established even for very 
Recently Smith et al. 


ligli bond strengths. 


63] combined a Gillespie type of approach with Green’s function 


reaction dynamics simulations to infer the effect of additional macromolecular crowders on 
virus capsid assembly. Here we have consider this important aspect in a fully spatial context. 
When placing non-reactive macromolecular crowding agents inside our simulation volume, 
we observe the same trends with regard to bond strength (optimal assembly at intermediate 
bond strength and crowding of our simulation volume at high and low bond strength). In the 
case of extra crowding, our simulations require the use of lower rates of capsomere influx in 
order to establish a continuous production of complete capsids, thus bringing our simulations 
closer to the physiological situation. 

Comparing the yield of complete capsids under static conditions with a fixed concentra¬ 
tion to the yield of complete capsids when gradually increasing the concentration, we demon¬ 
strated that the vulnerability of the assembly process to kinetic trapping can be significantly 
reduced if the concentration of capsomeres is dynamically increased. This conclusion agrees 
with the recent finding of a discrete stochastic simulation for genome-stabilized virus as¬ 
sembly that a linear increase in protein concentration dramatically increases the robustness 
against kinetic trapping [43]. 

To rationalize these results, it is helpful to think of assembly in terms of a free energy 
landscape (similar to transition networks in protein folding [6J, [65|). In the picture of 
the free energy landscape a kinetically trapped assembly process has reached a stable local 
minimum that prevents the formation of the desired minimum energy configuration. At high 
bond strengths and high concentrations the free energy landscape of the assembly process 
is very rough, and trapping in a local minimum is thus very likely. The gradual influx of 
new assembly material can be understood as a tilting of the free energy landscape. Thus, 
by continuously providing new capsomeres the assembly process can be guided towards the 
global minimum corresponding to the formation of full capsids and the dynamic capsomere 
supply prevents the system from becoming trapped in a local minimum configuration. 

As discussed before in our simulations we use enhanced assembly dynamics and relatively 
large influx rates to achieve reasonable computing times for our particle-based simulations of 
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empty capsid assembly. Thus our simulations include the full effect of diffusional encounters 
and excluded volume interactions. Our simulations suggest that qualitatively similar results 
are to be expected for lower rates of capsomere influx, which we use when accounting for the 
presence of additional macromolecular crowders. However, further progress in this direction 
needs algorithmic advances, including the use of GPU-code and analytical or resampling 
techniques to speed up simulation times 


66 


69|. Such advances then would allow us to 


also address more complicated virus architectures with different bond types, virus misfits, 
genome-assisted assembly and the interplay between virus assembly and gene expression. 

Although our simulations are performed with enhanced assembly dynamics and high rates 
of capsomeres influx, we believe that our findings have strong implications for the assembly 
process of virus shells under (dynamic) in vivo conditions. In particular our simulations 
suggest that kinetic trapping, which is often thought of as a major limitation preventing 
the successful assembly of complete capsids under (static) in vitro conditions, might only 
play a minor role in vivo if the supply of new capsomeres by the host cell is slow enough. 
Although assembly is still most efficient at intermediate bond strength, a dynamic supply of 
capsomeres should allow for robust self-assembly in a wide range of bond strengths without 
the need for additional helper proteins or scaffolds. 


results by Dykeman et al. 


'his finding is compatible with the 
431] and Hagan et al. 42] using different setups. Moreover, 


our dynamic, particle-based simulations show that self-crowding of the simulation volume 
can prevent steady-steady capsid assembly at high rates of capsomere influx and high bond 
strengths. 

Our simulations suggest that in the presence of other macromolecular crowders the effect 
of self-crowding due to the continuous production of capsomeres might occur for much lower 
rates of capsomere influx (and also lower capsomer concentrations). Here further investiga¬ 
tions are necessary to clarify the relevance of this regime for in vivo capsid assembly. One 
way to test our predictions is to use in vitro experiments. The setup studied in our sim¬ 
ulations could be experimentally realized using a microfluidic device that allows to control 
capsomere influx. At the same time, the bond strength might be controlled by changing tem¬ 
perature or ionic conditions. Capsid removal could be implemented simply by sedimentation 
42], by filtering or by boundaries that are sticky to completed capsids. For such a setup, 
we expect that the yield of full capsids at high bond strengths (low temperatures) depends 
on the rate of concentration increase. In particular, decreasing the rate of capsomere influx 
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should result in a higher yield of full capsids. 


As the assembly of a simple icosahedral capsid can be considered as a paradigm for protein 


assemblies 


70] or artificial assembly systems in general 


H 0, 


our findings also apply for 


other complex assembly structures. While recent advances in the design of artificial assembly 
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from colloidal particles have aimed at a dynamic control of the inter-particle interactions 
0 to increase the yield of the desired target structure, our simulations strongly suggest 
that a dynamic control of the material available for the assembly process can further help 
to increase the yield of the desired structure in artificial self-assembly systems. In these 
systems the effect of self-crowding discussed in our manuscript might also play an important 
role, depending on the rate of supply of new assembly material and the accessible assembly 
volume. 
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